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An efficient method is described to handle mesh indexes in multidimensional problems like 
numerical integration of partial differential equations, lattice model simulations, and determination 
f"^ ' of atomic neighbor lists. By creating an extended mesh, beyond the periodic unit cell, the stride 

in memory between equivalent pairs of mesh points is independent of their position within the cell. 
CNJ . This allows to contract the mesh indexes of all dimensions into a single index, avoiding modulo and 

5_j ' other implicit index operations. 
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PACS: 02.70.-c, 02.70. Bf, 05.10.-a, 45.10.-b 



Periodic boundary conditions are essential in all sorts of problems in solid state physics, condensed matter sim- 
ulations, and many other fields. The solution of partial differential equations by real-space discretizationEro, the 
calculation jof the interaction energy in lattice modelsa, and the determination of neighbor lists in molecular dynamics 
simulationso, are only a few of these problems. For illustration, let us consider the calculation of the laplacian of a 
^2 ■ function f(r) using finite differences. In three dimensions, one generally discretizes space in all three periodic direc- 
tions, using an index for each direction. For simplicity, let us consider an orthorhombic unit cell, with mesh steps 
Ax, Ay, Az. Then the simplest formula for the laplacian is 

^ fix,iy,iz (/<»+l,<y,«a ^fix,iy,iz fisc — l,iy,iz) I Ax 
+ (/ix^ + Mz ~2fi x ,i y ,iz + fix,i y -l,iz)l Ay 2 
~\~ {fi x ,i y ,i z + l ^fi x ,iy,iz fias , iy ,iz — 1 ) I Az 



A direct translation of this expression into Fortran90 code might read 



Lf(ix,iy,iz) = - f(ix,iy,iz) * (2/dx2+2/dy2+2/dz2) k 
+ ( f (modulo (ix+1 ,nx) , iy, iz) + f (modulo(ix-l ,nx) , iy, iz) ) / dx2 & 
+ ( f (ix, modulo (iy+1 ,ny) , iz) + f (ix,modulo(iy-l ,ny) , iz) ) / dy2 & 
+ ( f (ix, iy, modulo (iz+1 ,nz) ) + f (ix, iy ,modulo(iz-l ,nz) ) ) / dz2 

^ ' 

where the indexes i a (a — {x, y, z}) of the arrays f and Lf run from to n a — 1, as in C. There are two problems with 
this construction. First, the modulo operations are required to bring the indexes back to the allowed range [0, n a — 1]. 
And second, the use of three indexes to refer to a mesh point implies implicit arithmetic operations, generated by the 
compiler, to translate them into a single index giving its position in memory. 
' A straightforward solution to these inefficiencies would be to create a neighbor-point list j_neighb(i,neighb), of 
the size of the number of mesh points times the number of neighbor points. However, although the latter are only 
■ six in our simple example, they may frequently be as many as several hundred, what generally makes this approach 
unfeasible. A partial solution, addressing only the first problem, is to crcatc^ix (or more for longer ranges) one- 
dimensional tables Ja 1 {ia) — mod(i Q ±1, n a ) to avoid the modulo computations!]. Here, I describe a multidimensional 
generalization of this method, which solves both problems at the expense of a very reasonable amount of extra storage. 

The method is based on an extended mesh, which extends beyond the periodic unit cell, by as much as required 
to cover all the space that can be reached from the unit cell by the range of the interactions or the finite-difference 
operator. The extended mesh range is i™ m = —An a and i™ ax = n a — 1 + An a , where An a — 1 in our particular 
example, in which the laplacian formula extends just to first-neighbor mesh points. In principle, in cases with a small 
. £h ! unit cell and a long range, the mesh extension may be larger than the unit cell itself, extending over several neighbor 
cells. However, in the more relevant case of a large system, we may expect the extension region to be small compared 
to the unit cell. We then consider two combined indexes, one associated to the normal unit-cell mesh, and another 
one associated to the extended mesh 
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n ext _ jmax _ ^min _|_ _ n ^ _|_ 2An a . The key observation is that, if i ex t is a mesh point within the unit cell 
(0 < i a < n a — 1), and if j ex t is a neighbor mesh point (within its interaction range, i.e. \j a — i a \ < An a ), then the 
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arithmetic difference j ex t — i e xt depends only on the relative positions of i ex t and j ex t (i.e. on j a — i a ), but not on the 
position of i ex t within the unit cell. We can then create a list of neighbor strides Aij extl and two arrays to translate 
back and forth between i and i ex t- One of the arrays maps the unit cell points to the central region of the extended 
mesh, while the other one folds back the extended mesh points to their periodically equivalent points within the unit 
cell. Then, to access the neighbors of a point i, we a) translate i — > i ex t) b) find j ex t — i ex t + At jext] an( i c ) translate 
jext — ► j- Notice that several points of the extended mesh will map to the same point within the unit cell and that, 
in principle, a unit cell point j may be neighbor of i through different values of jext- In our example, the innermost 
loop would then read 

Lf(i) = 

do neighb = l,n_neighb 

j_ext = i_extended(i) + ij_delta(neighb) 

j = i_cell(j_ext) 

Lf (i) = Lf (i) + L(neighb) * f (j) 
end do 

where the number of neighbor points would be n_neighb=7, including the central point itself, and 
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Notice, however, that the above loop is completely general, for any linear operator, using an arbitrary number of 
neighbor points for its finite difference representation. In fact, it is even independent of the space dimensionality. 
Furthermore, the index operations required are just one addition and three memory calls to arrays of range oneH. 
This inner loop simplicity comes at the expense of the two extra arrays i_extended and i_cell (of the size of the 
normal and extended meshes, respectively) which are generally an acceptable memory overhead. Notice, however, 
that the the neighbor-point list ij_delta is independent of the mesh index i, what makes this array quite small in 
most problems of interest. 

Although particularly convenient for periodic boundary conditions, the present method may also be useful with 
fixed boundary conditions by, for example, mapping the outside points of the extended mesh to a single memory 
address with a fixed boundary value. In conclusion, I have presented an efficient method to handle index references 
in many typical problems involving periodic boundary conditions in more than one dimension. 
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